void Navier::NewFluidCellInitializer(int i, int j, int k)//Inicializador das c�lulas que eram Air e passaram aser fluido
{  
   //cout<<"\n Entrei na NewFluid \n";
   //cout<<"\n Init i,j,k ="<<i<<","<<j<<","<<k<<"\n";
   double tempmx, tempMx, tempmy, tempMy, tempmz, tempMz;
   int cardx, cardy, cardz, cardinal; 
   double   div = 0.0;
   unsigned int flag_level_temp = flag[i][j][k];
   tempmx=tempMx=tempmy=tempMy=tempmz=tempMz=0.0;
   cardx=cardy=cardz=cardinal=0;
   if(flag[i-1][j][k]>flag_level_temp) //&& (flag[i+1][j][k]<AIR_FLUID_END))
   {//Garantindo a divergencia nula: extrapolando a componente relativa a face "comum" e mantendo (repetindo o valor das outras)
      tempmx += U[i][j][k];//repetindo
      tempMx += 2*U[i][j][k]- U[i-1][j][k];//extrapolando
      tempmy += V[i-1][j][k];
      tempMy += V[i-1][j+1][k];
      tempmz += W[i-1][j][k];
      tempMz += W[i-1][j][k+1];
      cardx +=2;//indicador de c�lulas vizinhas por "face YZ" que s�o fluido; +2 indica que estou tratando a c�lula i-1,j,k
      cardinal++;//contabiliza a c�lulas vizinhas a c�lula i,j,k que s�o fluido
      double dif = tempMx + tempMy + tempMz - (tempmx + tempmy + tempmz)- div;
      /*if((fabs(dif) > 1E-6))
      {
         cout<<"\n\n ****divergente grande em (x,y,z) = ("<< i<<","<<j<<","<<k<<")";
         cout<< "\t flagi-"<< flag[i-1][j][k] << " dif = " << dif<<"\n";
      }*/
      div += dif;
   }
   //else if(flag[i-1][j][k]==flag_level_temp)
   else if(flag[i-1][j][k] >= AIR_FLUID_MIDDLE )
   {
     cardx +=2;
   }      
   if(flag[i+1][j][k]>flag_level_temp) //&& (flag[i-1][j][k]<AIR_FLUID_END))
   {
      tempmx += 2*U[i+1][j][k]- U[i+2][j][k];//
      tempMx += U[i+1][j][k];
      tempmy += V[i+1][j][k];
      tempMy += V[i+1][j+1][k];
      tempmz += W[i+1][j][k];
      tempMz += W[i+1][j][k+1];
      cardx++;// 1 indica que estou tratando a c�lula i+1,j,k; 3 indica que estou tratando as c�lulas i-1 e i+1
      cardinal++;
      double dif = tempMx + tempMy + tempMz - (tempmx + tempmy + tempmz)- div;
      /*if((fabs(dif) > 1E-6))
      {
         cout<<"\n\n ****divergente grande em (x,y,z) = (" << i << ","<< j << "," << k << ")";
         cout<< "flagi+"<< flag[i+1][j][k] << " dif = " << dif << "\n";
      }*/
      div +=dif;
   }
   //else if(flag[i+1][j][k]==flag_level_temp)
   else if(flag[i+1][j][k]>=AIR_FLUID_MIDDLE)
   {
     cardx ++;
   }      
   if(flag[i][j-1][k]>flag_level_temp) 
   {
      tempmx += U[i][j-1][k];
      tempMx += U[i+1][j-1][k];
      tempmy += V[i][j][k];
      tempMy += 2*V[i][j][k]- V[i][j-1][k];
      tempmz += W[i][j-1][k];
      tempMz += W[i][j-1][k+1];
      cardy +=2;
      cardinal++;
      double dif = tempMx + tempMy + tempMz - (tempmx + tempmy + tempmz)- div;
      /*if((fabs(dif) > 1E-6))
      {
         cout<<"\n\n ****divergente grande em (x,y,z) = (" << i << ","<< j << "," << k << ")";
         cout<< "flagj-"<< flag[i][j-1][k] << " dif = " << dif<<"\n";
      }*/
      div +=dif;
   }
   //else if(flag[i][j-1][k]==flag_level_temp)
   else if(flag[i][j-1][k]>=AIR_FLUID_MIDDLE)
   {
     cardy+=2;
   }
   if(flag[i][j+1][k]>flag_level_temp) 
   {
      tempmx += U[i][j+1][k];
      tempMx += U[i+1][j+1][k];
      tempmy += 2*V[i][j+1][k]- V[i][j+2][k];
      tempMy += V[i][j+1][k];
      tempmz += W[i][j+1][k];
      tempMz += W[i][j+1][k+1];
      cardy ++;
      cardinal++;
       double dif = tempMx + tempMy + tempMz - (tempmx + tempmy + tempmz)- div;
     /* if((fabs(dif) > 1E-6)||((i == 24) && (j==3) && (k==14)))
      {
         cout<<"\n\n ****divergente grande em (x,y,z) = (" << i << ","<< j << "," << k << ")";
         cout<< "flagj+"<< flag[i][j+1][k] << " dif = " << dif<<"\n";
      }*/
      div +=dif;
   }
   //else if(flag[i][j+1][k]==flag_level_temp)
   else if(flag[i][j+1][k]>=AIR_FLUID_MIDDLE)
   {
     cardy++;
   }
   if(flag[i][j][k-1]>flag_level_temp) //&& (flag[i][j][k+1]<AIR_FLUID_END))
   {
      tempmx += U[i][j][k-1];
      tempMx += U[i+1][j][k-1];
      tempmy += V[i][j][k-1];
      tempMy += V[i][j+1][k-1];
      tempmz += W[i][j][k];
      tempMz += 2*W[i][j][k]- W[i][j][k-1];
      cardz +=2;
      cardinal++;
       double dif = tempMx + tempMy + tempMz - (tempmx + tempmy + tempmz)- div;
      /*if((fabs(dif) > 1E-6))
      {
         cout<<"\n\n ****divergente grande em (x,y,z) = (" << i << ","<< j << "," << k << ")";
         cout<< "flagk-"<< flag[i][j][k-1] << " dif = " << dif<<"\n";
      }*/
      div +=dif;
   }
   //else if(flag[i][j][k-1]==flag_level_temp)
   else if(flag[i][j][k-1]>=AIR_FLUID_MIDDLE)
   {
     cardz+=2;
   }
   if(flag[i][j][k+1]>flag_level_temp) //&& (flag[i][j][k-1]<AIR_FLUID_END))
   {
      tempmx += U[i][j][k+1];
      tempMx += U[i+1][j][k+1];
      tempmy += V[i][j][k+1];
      tempMy += V[i][j+1][k+1];
      tempmz +=  2*W[i][j][k+1]- W[i][j][k+2];
      tempMz += W[i][j][k+1];
      cardz ++;
      cardinal++;
       //double dif = tempMx + tempMy + tempMz - (tempmx + tempmy + tempmz)- div;
      /*if((fabs(dif) > 1E-6))
      {
         cout<<"\n\n ****divergente grande em (x,y,z) = (" << i << ","<< j << "," << k << ")";
         cout<< "flagk+"<< flag[i][j][k+1] << " dif = " << dif<<"\n";
      }*/
   }
   //else if(flag[i][j][k+1]==flag_level_temp)
   else if(flag[i][j][k+1]>=AIR_FLUID_MIDDLE)
   {
     cardz++;
   }
   
   if(cardinal == 0)
   {      
      std::cout<<"no fluid neighbors-("<<i<<","<<j<<","<<k<<")\n";
      Initializer_No_Face_or_Edge_Neighbor( i, j, k);
      return;
      //exit(0);
   }

   if(cardinal == 6)
   {      
      std::cout<<"no freedom degree -("<<i<<","<<j<<","<<k<<")\n";
      flagsixfluidneighbors = true;
      return;
      //exit(0);
   }
   
   tempmx /=cardinal;
   tempMx /=cardinal;
   tempmy /=cardinal;
   tempMy /=cardinal;
   tempmz /=cardinal;
   tempMz /=cardinal;
   double deltax = tempMx - tempmx;
   double weightx = fabs(tempMx + tempmx);
   double deltay = tempMy - tempmy;
   double weighty = fabs(tempMy + tempmy);
   double deltaz = tempMz - tempmz;
   double weightz = fabs(tempMz + tempmz);
   double  excessx = 0.0; //N�o nulo apenas quando as duas faces "YZ" (paralelas) separam a c�lula vigente de c�lulas fluidos
   double  excessy = 0.0;
   double  excessz = 0.0;
   double  compensationx = 0.0; //Diferen�a a ser aplicada a componente u, devido ao excesso nas outras componentes
   double  compensationy = 0.0;
   double  compensationz = 0.0;
   
   //???
   double distributed_excess_factor_x=0;
   double distributed_excess_factor_y=0;
   double distributed_excess_factor_z=0;
   //fator relacionado a existencia de componentes, cujo valor nas duas faces
   //tranversais adjacentes a c�lula, j� foram computados. Se houver s� uma componente com essa propriedade o fator
   //vale 0.5 caso haja duas o fator vale 1, ou seja, com duas o valor � distribuido, com uma �nica o valor � compensado
   //em �nica dire��o.
   if(cardx ==3) 
   {
      excessx = deltax - U[i+1][j][k] + U[i][j][k];
      if(cardy ==3)
      {
         excessy = deltay - V[i][j+1][k] + V[i][j][k];
         distributed_excess_factor_z =1.0;
      }
      else
      {
         if(cardz ==3)
         {
          excessz = deltaz - W[i][j][k+1] + W[i][j][k];
          distributed_excess_factor_y =1.0;
         }
         else
         {
            double sum_weight = weightz+weighty;
            if( sum_weight < 1E-8)
            {
               distributed_excess_factor_z = 0.5; 
               distributed_excess_factor_y = 0.5;
            }
            else
            {
               distributed_excess_factor_z =weightz/sum_weight; 
               distributed_excess_factor_y =weighty/sum_weight; 
            }
         }
      }
   }
   else
   {
      if(cardy ==3) 
      {
        excessy = deltay - V[i][j+1][k] + V[i][j][k];
        if(cardz ==3)
        {
            excessz = deltaz - W[i][j][k+1] + W[i][j][k];
            distributed_excess_factor_x = 1.0;
        }
        else
        {
            double sum_weight = weightz+weightx;
             if( sum_weight < 1E-8)
            {
               distributed_excess_factor_z = 0.5; 
               distributed_excess_factor_x = 0.5;
            }
            else
            {
               distributed_excess_factor_z =weightz/sum_weight; 
               distributed_excess_factor_x =weightx/sum_weight; 
            }
        }
     }
     else
     {
        if(cardz ==3)
        {   
            excessz = deltaz - W[i][j][k+1] + W[i][j][k];  
            double sum_weight = weighty+weightx;
             if( sum_weight < 1E-8)
            {
               distributed_excess_factor_y = 0.5; 
               distributed_excess_factor_x = 0.5;
            }
            else
            {
            distributed_excess_factor_y =weighty/sum_weight; 
            distributed_excess_factor_x =weightx/sum_weight;
            }      
        }
     }
   }
   if(cardx <=2)
   {
      compensationx = (excessy +excessz) * distributed_excess_factor_x;
      if(cardx == 0)
      {
         U[i][j][k] = tempmx - compensationx * 0.5;
         U[i+1][j][k] = tempMx + compensationx * 0.5;
      }
      else
      {
         compensationx +=deltax;
         if(cardx == 1)
            U[i][j][k] = U[i+1][j][k] - compensationx;
         else
            U[i+1][j][k] = U[i][j][k] + compensationx;
      }
   }
   if(cardy <=2)
   {
      compensationy = (excessx +excessz) * distributed_excess_factor_y;   
      if(cardy == 0)
      {
         V[i][j][k] = tempmy - compensationy * 0.5;
         V[i][j+1][k] = tempMy + compensationy * 0.5;
      }
      else
      {
         compensationy +=deltay;
         if(cardy == 1)
            V[i][j][k] = V[i][j+1][k] - compensationy;
         else
         V[i][j+1][k] = V[i][j][k] + compensationy;
      } 
   }
   if(cardz <=2)
   {
      compensationz = (excessy +excessx) * distributed_excess_factor_z;   
      if(cardz == 0)
      {
         W[i][j][k] = tempmz - compensationz * 0.5;
         W[i][j][k+1] = tempMz + compensationz * 0.5;
      }
      else
      {
         compensationz +=deltaz;
         if(cardz == 1)
            W[i][j][k] = W[i][j][k+1] - compensationz;
         else
            W[i][j][k+1] = W[i][j][k] + compensationz;
      }
   } 
   //double div_1 = U[i+1][j][k] + V[i][j+1][k] + W[i][j][k+1] - (U[i][j][k] + V[i][j][k] + W[i][j][k]);
   
   /*if (div_1 > 1.0e-8)
   cout << "\n div 1: " << div_1;*/
   /*if(((j==1) && (k==7) && (i==25))||
      ((j==1) && (k==7) && (i==24))||
      ((j==1) && (k==8) && (i==25)))
   //if((cardx==3)|| (cardy==3) || (cardz==3))
                 {
                    cout<<"\n i,j,k ="<<i<<","<<j<<","<<k<<"\n";
                    cout<<"cardinal,cardx,cardy,cardz ="<<cardinal<<","<<cardx<<","<<cardy<<","<<cardz<<"\n";
                    cout<<"tempmx,tempmy,tempmz ="<<tempmx<<","<<tempmy<<","<<tempmz<<"\n";
                    cout<<"tempMx,tempMy,tempMz ="<<tempMx<<","<<tempMy<<","<<tempMz<<"\n";
                    cout<<"\n U,V,W new ="<< U[i][j][k]<<","<<V[i][j][k]<<","<<W[i][j][k]<<"\n";
                    cout<<"\n U,V,W new ="<< U[i+1][j][k]<<","<<V[i][j+1][k]<<","<<W[i][j][k+1]<<"\n";
                 }
   //cout<<"\n Sai da NewFluid \n";*/
}